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ABSTRACT 

The Reuven Ramaty High Energy Spectroscopic Solar Imager (RHESSI) uses the rotational modulation 
principle (Schnopper et al., 1968) to observe temporally, spatially, and spectrally resolved hard X ray and 
gamma ray images of solar flares. In order to track the flare evolution on time scales that are commensurate 
with modulation, the observed count rates must be demodulated at the expense of spatial information. The 
present paper describes improvements of an earlier demodulation algorithm, which decomposes the observed 
light curves into intrinsic source variability and instrumental modulation. 



INTRODUCTION 

The RHESSI (Lin et al. 2003) instrument is designed to explore particle acceleration in solar flares (e.g.. 
Miller et al., 1997). It observes photons at energies from 3 keV to 17 MeV with spectral resolution up 
to £'/Aii^~500. In order to obtain images, RHESSI is equipped with 9 pairs of X ray shadowing grids (9 
'subcollimators'), which are fixed on the rotating spacecraft (spin period rs'~4s). The instantaneous trans- 
mission probability of one subcollimator, projected onto the solar disc (heliocentric cartesian coordinates 
x), is called the 'modulation pattern' (Hurford et al., 2003a). It has the approximate form (i=1...9) 

M^x, t) = 4 + ai cos {k\t) ■ (x - P(t)) + V') . (1) 

The wave vectors k*(t) correspond to periods ('angular pitches') of ~2.6 • 3*/'^ asec, and P(t) is the imag- 
ing (optical) axis. The angle ^{^{t) between k*(t) and the heliocentric x-axis is given by ^]^{t) = 7r/2 — 
^roii(i) — '^grid' where <I>roii(*) is the roll angle and $gj.id is the grid orientation on the spacecraft. The triples 
(<I>roii(i)) P(i)) are referred to as 'aspect solution', and are continuously monitored by the on-board aspect 
systems (Fivian and Zehnder, 2003; Hurford et al., 2003b). The coefficients aQ,a\,'ip^ in eq. (Q) depend on 
energy and (weakly) on x — P(t); they describe the X ray transmission of the grids. Above 100 keV the grids 
become increasingly transparent. In what follows, the energy band is fixed, and the energy dependence is 
omitted. 

At any time, the 9 subcollimators respond to 9 different Fourier components of the solar brightness 
distribution. During half a RHESSI rotation, up to thousand Fourier components are collected, from which 
the spatial brightness distribution can be restored (Hurford et al., 2003a). The imaging information is thus 
encoded in the temporal modulation of the observed count rates. In order to track impulsive events (energy 
release, acceleration) on time scales below 2s, the observed count rates must be demodulated at the expense 
of spatial information. 

The present article reports on further development of the demodulation method of Arzner (2003), which 
compared count rates with similar modulation patterns. This constraint was rather restrictive, since both 
the roll angle and the grid phase had to agree. The request for grid phase agreement can in fact be relaxed 
by the use of visibilities^; this is described in Section 2 below. In order to work properly, this approach 

^At each time, the two 'visibilities' are the projections of the instantaneous brightness distribution onto sine- and cosine 
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Fig. 1. Left: Motion of the imaging axis P(t) (circles) during Np = 8 spacecraft revolutions at the Flare of 
August 23, 2002, 06:25 UT. The wave vectors (arrow indicates orientation of at final time) rigidly co-rotate 
with P(t). The average spin axis is marked by a cross. Right: change of modulation phase of subcollimator 3 
per At at the average spin axis {Acj) = cpAt). The condition \A(j)\ <^ In is the basis of Eq. (jH)). 

requires that the time binning is synchronized with the spacecraft rotation, i.e., that <&roii (mod vr) recurs 
after an integer number of time bins At. On the other hand, At must be an integer multiple of a binary 
micro second, which is the hardware time resolution. The improved version choses an optimum compromise 
between the two conflicting constraints, while keeping the number of counts per bin at a statistically uselful 
level (~ 10). Section 3 shows benchmark simulations and the demodulation of real RHESSI data. 

DEMODULATION METHOD 

The goal of the method is to separate spatial and temporal variations of the solar brightness distribution. 
The basic observation is that the visibilities at equal grid orientation should be compatible if the brightness 
distribution did not change with time. Deviations are attributed to time dependence of the source distri- 
bution. Technically, this is done by a regularized maximum-likelihood fit with many (>10^) parameters 
describing both spatial and temporal source variability. 

Let us start with a brief description of the spacecraft motion, which determines the time-dependent 
instrumental response (Eq. d). During a few {Np < 10) spin periods, the RHESSI motion is composed of 
an approximately uniform rotation and a slow translation (Figure ^ left) . We formalize this by making the 
following assumptions (which are checked from the aspect solution): (i) the spin period Ts is constant to 
accuracy At over the time of interest (NpTs); (ii) the wave vectors k*(t) rotate clockwise with spin period 
Ts] (iii) the imaging axis P(t), which is generally not aligned with the spin axis^ s(t), rotates (clockwise, 
period Ts) around s(t). (iv) the spin axis and |s(t)-P(t)| vary slowly due to precession {fprec ~ 0.015 Hz), 
inertia adjustment, and magnetic torquing (Lin et al., 2003). 

We next set up a model for the solar brightness distribution. It is assumed that the true brightness can 
be represented as 

5(x,t) =^rfc(t)5fc(x) (2) 
k 

where -Bfc(x) are arbitrary non-negative normalized (/ 5fc(x)(fx=l) functions and rfc(t) (ct/s) is the instan- 
taneous count rate. In principle, Eq. (jSJ describes all possible brightness distributions, but in practice we 
restrict ourselves to a few source components k; perhaps, a gradual plus an impulsive component with dif- 
ferent intrinsic time scales. Impulsive components could be associated with magnetic footpoints. It should 
be pointed out that Eq. @ cannot account for continuously moving sources, and that non-solar background 
is neglected. 

Combining the RHESSI response (Eq. ^ with the brightness model (Eq. [2} one can predict the expected 
counts in time bins At. As we wish to resolve the source evolution it is assumed that rfc(t) varies slowly 

components of the modulation pattern. 

^Strictly speaking, s{t) denotes the penetration point of the instantaneous RHESSI spin axis through the solar disk. 



3 



during At. The expected number of counts in the i-th subcolhmator and t-th time bin is then given by 

\l = LlJ2rk,tml, (3) 

k=0 

where the hvetime measure LI accounts for detector dead time (Smith et al., 2003) and data gaps^, r^^t = 
(At)^^ J^^ ^ dt' rk{t') is the (discrete) count rate, and 

mU= / dt' dxM\-K,t')Bk{-K). (4) 



We recall that x denotes heliocentric coordinates. While being natural, these coordinates are not well 
adapted to the observing geometry, and do not allow to separate fast (rotational) and slow (translational) 
contributions to modulation. Such a separation is, however, possible if we set x = x' + (s), where (s) is the 
average spin axis (Figure Q left cross; average over Np spin periods). After expanding the cosine in Eq. 
Eq. ((3} becomes 



k,t 



^ dt' j d-x! (^4{t') + ai(t') cos(k^(t') • x') cos {k^t') • ((s) - P(t')) + V''(i'))} 

-al(t') sin(k*(t') • x') sin {k^(t') • ((s) - V(t')) + V'^(t')}) i?fe(x' + (s)) . (5) 



Now, the term in curly brackets is the grid phase at the average spin axis, which is a slowly varying quantity. 
If also ao(t) and a\{t) are weakly varying, then Eq. © can be approximated by 



/•t+At /• 

ml ^ ~ a\^^^At + a\^^ cos </.J / dt' / dx' cos(k*(t') • x')Bk{x' + (s)) 
/•t+At /• 

J dt' J dx! sin(k*(t') • x')5fc(x' + (s)) (6) 



—"i.t ^^^^ 

<At + a\ t cos CU - a\, sin 0J SU (7) 



where ^ is a discrete version of a\{t), (j)\ = \i\ - ((s) — Y't) + ipl is the discrete grid phase at average spin axis, 
and the last line defines the (At-integrated) visibilities ^ and ^. The visibilities satisfy (C|, ^■)^ + (5| j)^ 
< (At)^, ~ ^fcj' ^kj+Ns/2 ~ ~^kj-> where Ns is the (even) number of time bins per spin 

period. The visibility approach (Eqns. El and 1^1 is valid if \rk/rk\At <C 1 and 



|d*o/aol^* < 1' \d\/a\\At < 1, (27r)-^|(/)*|At < 1 , (8) 

which must hold during the whole time interval NpTg. A safe threshold for Eqns. (jSJ is 0.1, which is checked 
by the computer code. This is usually uncritical (Figure fright): from typical aspect data, one finds that 
time bins as large as At = 0.3s are admissible even for the finest subcolhmator. It should be stressed that 
this is much larger than the modulation period - the modulation must not be resolved to predict similar 
(At-integrated) visibilities. Of course, integrating over modulation razes spatial information, but this does 
not affect the time profiles r^^f 

At this point we have completed the forward model (Eqns. I3I7|1 . It has parameters {C|, j-, S], ^j, where 
1 < i < 1 < t < Ns-Np, and i is out of a subcolhmator set satisfying Eq. 0. We turn now to 

the counting statistics. Let the observed light curve of subcolllimator i have cj counts in time bin t. These 
counts should scatter around the expectation value (Eq. according to Poisson statistics. The agreement 
between observation and prediction is therefore measured by the Poisson log likelihood ratio (Eadie et al., 
1971) ('C-statistic') 

logL = 5:*{-Aj + cj(l + ln4)}. (9) 



it 



-t 



^Data gaps occur independently in different subcollimators and have durations of milliseconds to seconds (Figure]^). They 
affect some 30% of all data, and are presumably caused by cosmic rays. The livetime measures L] are between zero and one. 
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The likelihood L is proportional to the probability that {cj} is observed if {A^} was true. Log L is negative, 
and reduces to —^X^ iii the limit of high count rates. The asterisks in Eq. Q indicates that the sum is 
restricted to count rates with livetimes LI > L^i^, with Lmin typically chosen as 0.5. This prevents from pile 
up and from redundant numerical operations, but Li = is frequent ( > 25% 

data gaps). The rejection of low livetimes is not expected to affect the result. 

Since there are more model parameters {CI j, j,rk^t} than observations {cj}, and since the model 
parameters are not independent of each other, the maximum-likelihood problem does not have a unique 
solution. In order to regularize the problem we assign an a priori probability exp{— ^ J2ti^k,t+i — '^fc,*)^} to 
the k-th time profile r^^f This choice favors smooth time profiles, and allows an explicit trade-off between the 
smoothness of the demodulation and the likelihood. It also has the useful feature that it linearly interpolates 
across data gaps. Taking the a priori probability into account, the logarithm of the total (observational -|- 
a priori) probability P becomes, according to Bayes' rule, 

log P = log L - i ^ ak{rk,t - rk,t-if ■ (10) 

kt 

The quantity log P possesses a unique maximum which is found by Newton/Marquardt type iterations. The 
coefficients determine the smoothness of the solutions rt^k- Small Ok yield small correlation times of 
rt^k- While the exact relation between and r/^ depends on the actual data set and on the definition of r^, 
there exists a useful empirical estimate, 

Tk ^ max{At,^/a^{c)/{aoL)) . (11) 

In Eq. (|11() . Tfc is defined via the normalized autocorrelation Ak{t) of r^^fc by Aj.{t) ~ 1 — ^{t/rj^)^ + 
O(i^). The accuracy of Eq. (|11|) is within a factor two. The adjustment of ak is done by the user. The 
Newton/Marquardt iterations have the property that two solutions rt^k agree if their ak agree. 

RESULTS AND DISCUSSION 

As a benchmark, Figure|2K-e shows simulations where the model assumption (Eq. [2} is met. There are two 
sources at (-715", 630") and (-650", 690"), of size 1" and 3", with impulsive (Figure[2^ solid line) and gradual 
(Figure [2^ dashed) time profiles. The small source sizes yield maximal modulation and therefore represent 
a worst-case test for demodulation and residual modulation. Fig[2l3 shows a close-up of the true impulsive 
time profile, which is to be compared with the estimates c)-e) obtained from simulated 'observations' with 
{cl)/At = 6 • 10^. The simulated aspect solution has the form P(t) = Vpreci + ?^cone(cos((5 — Ot), sin((5 — Qt)) 
with Vprec = (0.8,1.3) asec/s the instantaneous precession rate, = 27r/4s the RHESSI angular frequency, 
Tcone = 120" the couiug radius, and 6 a phase offset. Figure HJ^-d shows the demodulation results (solid line) 
and averages over subcollimators 1-3 and time (dotted line) with similar nominal resolution. In Figure [2t) 
At=5ms, and the demodulation ro^t+Ti,t (solid line) of subcollimators 1-9 is compared with a moving average 
{c)/{aQL) (dotted) over subcollimators 1-2 and time interval Ta=0.1s. The autocorrelation time constant of 
the impulsive component ri^t is ti = 0.09s, while Eq. (|11() gives ri ~ 0.14s. Numerically, the rms deviations 
from the true time profile are 5% (average) and 3% (demodulation), while the statistical error assigned to the 
At-binned true light curve would be about 6%. While the benefit from 3% to 5% is rather small and within 
the statistical error, the situation changes if data gaps are present. Figure |2jl is similar to [St, but includes 
simulated data gaps with {L\)=Q.7. Simultaneous data gaps in subcollimators 1-3 with duration >Ta cannot 
be interpolated with a moving average (dotted line), and the rms deviation of the moving average from the 
true time profile is 25% (whereas it is 4% for the demodulation). If one tries to improve the moving average 
by including coarser subcollimators one encounters the difficulty that their (at least, slowest) modulation 
commensurates with the time scale to be resolved. The moving average then fails even in the absence of data 
gaps (Figure 12^ dotted). A possible way out is demodulation as described in this paper (Figure [2^ solid). 
Here, both the demodulation (solid) and average (dotted) involve subcollimators 7-9. The rms deviation 
from the true time profile is 6% (demodulation) and 18% (average), respectively. A test was also made 
for violation of the model assumption ^ . Figure |2f shows the demodulation of (simulated) counts from 5 
sources with different time profiles. The demodulation (Figure [3 middle curve) uses subcollimators 1-9 and 
is based on the (wrong) assumption that -B(x, t)=ro(t)i?o(x)+ri(t)i?i(x). Nevertheless, it has slightly better 




Fig. 2. Benchmark simulations: a) true time profiles of impulsive (solid) and gradual (dashed) sources; b) close- 
up; c) demodulation from subcollimators 1-9 (solid), compared to an average over subcollimators 1-3 and time 
Ta=0.1s (dotted, shifted for better clarity); d) similar to c) including 30% simulated data gaps; e) demodulation 
(solid) and average (dotted, ra=0.1s) of subcollimators 7-9 (no data gaps), f): robustness against violation of 
Eq. (121) - top to bottom curves: true unmodulated light curve of 5 independent components; two-component 
demodulation from subcollimators 1-9; average over subcollimators 1-3 and time ra=0.25s. The middle and 
bottom curves are shifted for better clarity, g): iterative convergence of Eq. ()10jl: — log P (boldface); — log L 
(solid); gradual (^=0, dotted) and impulsive {k=l, dashed) smoothness constraints. 

rms deviation from the true (spatially integrated) light curve (Figure [2f top) than a corresponding average 
(Figure[2t dotted line). Surprisingly, violation of the model assumption (jSJ does not degrade the result below 
the quality of a simple average. The iterative convergence to the maximum-P solution is demonstrated in 
Figure Ef, referring to simulation [2t) . The total probability (boldface) is dominated by the log likelihood 
(solid line), indicating that the demodulation is primary determined by the agreement with the observation. 

The algorithm is then applied to RHESSI observations of solare flares. A flare with rich temporal fine 
structures occurred on August 23, 2002, 06:25 UT (Figure 01) • The energy band under consideration is 6-25 
keV, and the mean observed count rate is 1700/s/subcollimator. The time bins are 9.78ms, so that there are 
^ 10 counts per time bin. The flare center is at (920", -160") (Figure |211-e). The time scale achievable in the 
demodulation is in the order of 0.1s, by the order-of-magnitude argument that for each (temporal -l-spatial) 
degree of freedom there should be some 100 photons. A first estimate on the true light curve is provided by 
the average {c)/{aoL) over subcollimators 1,3,4 and time ra=0.08s (Figure |3K-b, dotted, shifted for better 
clarity). The demodulation (solid line) is obtained from all subcollimators except #2 and #5, which are 
excluded because of high background (#2) (Smith et al., 2003), and a distinctly inferior likelihood (#5) 
which is not fully understood at present. The decomposition into impulsive and gradual components (Eq. 
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Fig. 3. Flare of August 23, 2002, 06:25 UT: a) average over subcollimators 1,3,4 and time Ta=0.08s (dotted, 
shifted for better clarity), demodulation rQ^t+fi^t (solid line), and gradual component ro,f (smooth curve) from 
subcollimators 1,3,4,6,7,8,9; b) close-up of demodulation (solid) with errors (gray, see text) and moving average 
(dotted); c) example of observed (crosses) and predicted (gray line) counts for subcollimator 6. e), d): CLEAN 
images at total, gradual maxima (axes in heliocentric arc seconds). 



121) removes data gaps, and also shows that the gradual component appears delayed (Figure smooth 
curve) . This seems consistent with a spatial change of the brightness distribution (Figure Otl-e) , suggesting 
the emergence of a gradual (possibly thermal) source at ~(924" ,-155"). Figure Ot shows an example of 
observed counts (crosses), together with the predicted Poisson parameter (gray line). The data gaps have 
zero counts. The correlation time of ri^t derived from the autocorrelation is ri = 0.163s; the estimate (Eq. 
[TT|) gives 0.164s. An important issue is the reliability of the demodulated fine structures. This can be 
tested by locally perturbing the maximum-P solution ri j until log L changes by a given amount (Eadie et 
al. 1971). Adopting the conventional threshold AlogL = ^, one obtains the error bars shown in Figure 
1313. The average error is 1400 ct/s. Large error bars are characteristic for times with data gaps in several 
subcollimators. 

Let us finish this discussion with a last but important point. Although the formalism of Eq. H10() allows 
for arbitrarily small count rates, there is no gain in temporal resolution if At is made so small that only 
few (or fractions of) photons are contained in a time bin. At such low count rates, the Poisson error is 
large, and Eq. (\lU\i becomes dominated by the smoothness constraint. Higher time resolution can only be 
achieved at higher count rates, or at the expense of statistical significance {xied < As a conservative 
rule-of-thumb, a few 10^ counts/s are needed to resolve 100 ms structures on the 5% level. It is not easy 
to beat this limit, not even for extremely transient and bright events, because their modulation is hard to 
estimate and dominates the Poisson error. 
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